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SUMMARY 


A method for analyzing viscous nonadiabatic flow inside 
turbomachine passages has been presented earlier in Vol. I of 
this study. The reported analysis dealt with incompressible 
laminar flow inside straight radial bladed impellers. In the 
present volume the developed field analysis is expanded to study 
turbulent flow within turbomachines having arbitrary blade 
geometries. Effects of turbulence are modeled using two equations, 
one expressing the development of the turbulence kinetic 
energy and the other its dissipation rate. To account for 
complicated blade geometries, the flow equations are formulated 
in terms of a nonorthogonal boundary fitted coordinate system. 

The analysis is applied to a radial inflow turbine. The 
solution obtained indicates the severity of the complex inter- 
action mechanism that occurs between the different flow regimes 
(i.e., boundary layers, recirculating eddies, separation zones, 
etc.). Comparison with nonvisceus flow solutions tend to justify 
strongly the inadequacy of using the latter with standard boundary 
layer techniques to obtain viscous flow details within turbo- 
machine rotors. Capabilities and limitations of the present 
method of analysis are discussed. 

The computer used in this work is an AMDAHL 470. The flow 
domain has been divided into 30 x 40 step sizes. Typical CPU 
time for the case studied is 520 seconds. 

INTRODUCTION 

During the last decade extensive progress has been made 
towards the development of suitable solution methods to the flow 
in turbomachines. Most of these methods dealt mainly with the 
flow field inside turbomachines as being inviscid [1, 2]. Few 
studies, on the ether hand, considered viscous effects using 
some kind of flow approximation 13, 4, 5). The general flow 
solutions within turbomachine passages, however, remains mostly 
beyond the scope of the prevailing modes. 

These limitations, especially of the numerical methods, 
have long been known to be acute due to the complicated nature 
of the flow' problem. The blade rows impose vorticity and 
velocity fields on the flow. The inviscid effects are dominant, 
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yet the viscous effects are not small enough to be neglected. 
Because such flows can only be described using the fully 
viscous equations of motion, the limitations of the prevailing 
methods is a direct consequency of at least three major 
obstacles. The first is the excessive computer time required 
for the solution of the equation of motion. The second is the 
difficulties and uncertainties associated with the treatment 
of certain boundary conditions. The third is the lack cf 
accurate and general methods for describing the turbulent 
transport phenomena that takes place inside tha blade rows. 

In an earlier report [7] , the authors of th:.s article 
reported a method suitable for analyzing the viscous incom- 
pressible flow in turbcmachines that partially overcome the 
first two obstacles mentioned above. In principle, the approach 
emphasized the feasibility of obtaining viscous flow details 
within turbomachine passages by appropriately combining several 
blade-to-blade viscous flow solutions. Each of these solutions 
were obtained through the numerical integration of the full 
Navier-Stokes equations over a predetermined computational 
surface that extends between the blades. The set of computa- 
tional surfaces required for the analysis (see Fig. 1) 
were themselves generated from the solution of the 
Euler ecruations in the nanner sugaested bv w u [l]. 

The results that had been obtained by applying this method of 
analysis to flows within straight radial bladed impellers were 
prc..iising. This was evidenced by the fair agreement achieved 
between the predicted and experimental data [7]. 

Recently, the above analysis was extended to deal with 
nonadiabatic incompressible viscous flows inside turbomachines 
having arbitrary blade geometries. In particular, the approach 
which has been established is modified and utilized in con- 
junction with: 

1. A two equation turbulence model. 

2. A coordinate transformation package [6] that employs 
a nonorthogonal boundary-fitted coordinate system to 
suit the most complicated blade geometries. 
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The differential equations that describe the viscous flow 
in turbomachine passages are briefly reviewed in this article 
with pertinent references made to the questions of their solution 
efficiency and accuracy of describing the different boundary 
conditions. The numerical procedures developed on this basis 
are discussed and illustrated by results obtained for the flow 
in a cooled radial inflow turbine. 

1. MATHEMATICAL FORMULATION 
Governing Differential F-~ uaticns 

The development of the mathematical model for the present 
problem starts with the Reynolds equations for the mean turbu- 
lent motion of flow. These equations are written in a generalized 
form for a turbomachine rotor as follows [7] : 

Conservation of Mass ; 

V- (pW) = 0 (1) 

Conservation of Momentum : 

D (W • VW + 2Q x W - P. 2 R) 

= - 7p - Vx[u e (VxW)] + j 7(u e V-W) (2) 

Conservation of Energy : 

0<W-Vh) = + ( W • 7 ) p + D + 7 • ( K7T) (3) 

where p, o, T, h and D denote the static pressure, 
density, temperature, enthalpy and dissipation function 
respectively. While, W represents the relative velocity vector 
at any point whose location is def'ned by the relative position 
vector R in a rotating frame of reference (Fig. 1). The speed 
of rotation of this frame is equal to the angular velocity of 
the machine (?.) . The effective turbulent viscosity in the 
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above equations is assumed to describe the effects of Reynolds 
stresses and is determined, in the present study, from a two 
equation turbulence model. 

As pointed out in the introduction and elsewhere [7] , the 
flow equations (1) to (3) are solved on several blade-to-blade 
computational surfaces. The geometry of such surfaces can be 
conveniently described in terms of the coordinate (m,$,n) as 
shown in Fig. 1. To reduce the complexity of handling the 
computations along these surfaces, it is found more convenient [7] 
to introduce a stream function and the mean vorticity of the 
turbulent flow as new dependent variables. The flow governing 
equations (1) through (3) , when recasted in terms of these 
variables and transformed to the (m,0,n) coordinate system result 
in the following system of equations [7] . 


Stream Function Equation : 

_ 1 r 9_ ,Mr 3<K 3_ , M jhh , 

r 3m op 3m 3<J> ^bpr 3 <tr 


(4) 


In this equation u denotes the mean vorticity and is the 
stream function defined according to the following relations: 


W 

m 


M 1 dji 
b pr 3$ 


and 


W 


M 1 3 i|i 
b p 3m 


(4a) 


where W , W. are the components of the mean relative velocitv 
m <p . 

vector W in the meridional, tangential directions, and M is 
the mass flow passing through the volume bounded by the surface 
and the filament thickness, b, shown in Fig. 1. The defi- 
nition of as given by equation (4a) identically satisfies 
the continuity equation (1) . In the process of eliminating 
the pressure by taking the curl of the momentum equation (2) , 
while using equation (4) and (4a), the following equation is 
obtained . 
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Vorticity Transport Equation ; 


3_ ,M 3'Ji ^ _ 3_ 
dm b 30 UJ ' 3$ 


M SO 
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+ G, 
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( 5 ) 


The source term in the above equation represents the generation 
or decay of vorticity due to the effects of rotation, it is 
given by; 


= 2CH 
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Energy Equation ; 

It is found appropriate for a turbomachine rotor to express 
the energy equation in terms of the total enthalpy (H) of the 
flow. The total enthalpy for a turbulent flow is expressed 
as follows: 

W 2 r. 2 r 2 

H = h + j- + ftW^r + — — + E (6) 

where E is the kinetic energy of turbulence. Thus, the energy 
equation (3) when transformed to the (m,0',n) coordinate system, 
results in [71 : 

M ,3_,„ 3jK _ £_/ h 3’K -I _ r 3 Hv _ 1_ _3_ ,^_e 3H 

b l 3m v 3<p 30' 3m 3m'?r 3nr r 3<p v Pr 3o ; 
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{u e rt 5F 

3W 2 /2 


‘s - 

CE 
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3m 

3m 
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l r L Pr 

3W 2 /2 


‘s - 

b CE 
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3 0 


Pr' 
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I— (u u) 
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dtp 

('•V 1 

- Dr 

4- G 


5 



where the source term G£ is given by 
G 


2 f<V + r2a) t ] -l? 1( V + r2n) Is 1 


+ fl{ i r Is + ^ 1 


13 , e 3 


fir' 


+ 7 tT Ip? fr (W.r + ^)] } 
r o<p Pr dip ip 2 


(7a] 


and the dissipation function# D, is expressed as 
3W 2 . 3W. W _ 

D * 2 « a {< 3iT> + <? + ST sin “> > 


3W, , 3W W. 0 

, r cp JL m ^ i ^ 

+ y i*^r— + — -r-— - — sina) 
e cm r 3b r 


(7b) 


In the above equations Pr denotes the turbulent Prandtl number 
and S„ is the turbulent Schmidt number for the kinetic energy 
of turbulence E. 


Turbulence Model 

The mean flow equations introduced earlier as equations (4) , 
(5) and (7) involve more unknowns than appropriate to the 
equations. This well known problem of closure is a character- 
istic of all nonlinear stochastic systems. The majority of 
investigators in the area of turbulent flow have been concerned 
with this closure problem. The most prominent and practical 
attempts have been the development of various second orde.’- 
two equation turbulent models. In this, two scalar quantities 
are used to characterize the turbulence, with each quantity 
satisfying a nonlinear partial differential equation. Commonly, 
the two equation models employ the kinetic energy of turbulence 
as one of the quantities. The second quantity used by different 
researchers [0] is a dissipation function, a length scale, or 
a dissipation rate. 
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The two equation model proposed by Spalding [9] for the 
kinetic energy of turbulence, E, and its volumetric rate 
of dissipation e is used in the present work. Previous 
experience [10] indicates that this model provides accurate 
prediction of various flow regimes within turbomachine passages, 
in addition to being less complex and more general. The 
complete details of the model and its implications are given 
in reference [9]; in brief, its main features are outlined 
in what follows. 

The Reynolds stresses are related to the velocity gradients 
through an effective viscosity, u . This viscosity is assumed 
to be isotropic at any point, but can vary throughout the flow 
field. The distribution of the effective viscosity u e is 
determined from the local values of density, o, the turbulent 
kinetic energy, E, and its volumetric dissipation rate, c. 

The expression used is: 


u 

e 


+ C D 



(7c) 


where the laminar viscosity, , is considered to be uniform 
and known . 

The equations governing the transport of E and e have 
a general form which includes convection, production, diffusion 
and dissipation terms. These can be written in terms of the 
(m,0,n) system of coordinates as follows [11]. 

Turbulent Kinetic Energy Equation : 


M 

E 



fp 2l2.\ _ 2 /p iLii\ 1 

(E 30 (E 3m 


3_ ( ^e_ 3E 

30 r30 


3 , J e aE. 
T— (r — r -p— ) 
om am 


) - rD + oer = 3 


( 8 ) 



D i • •" ipation Rate Equation : 
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b 
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+ C 


2 2 


pr = 0 


(9) 


The terms containing the coefficients S cc and S c _ in equations 
(8) and (9) represent the diffusion rates of E and c 
respectively. These coefficients therefore have the signifi- 
cance of turbulent Schmidt numbers for the transport process 
in question. The model contains five empirical constants 
which are assigned the values given in the following table. 


Values of the Empirical Constants 
for the k-e Model of Turbulence 



0.09 1.44 1.92 1.0 1.3 


The basis for choosing the above cited constants are discussed 
in details in reference [9]. It may suffice therefore to remark 
that C, is chosen so that Von Karman constant equals 0.42 and 
C, is determined by reference to the decay of grid turbulence. 

C D is fixed from the following requirements in a constant stress 
layer : 



where t is the wall shear stress. The diffusion coefficients 

V* 

S rT , and S c were fixed in [9] by computer optimization. 
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It should be emphasized at this point that the additional 
production cf turbulence energy by the rotational effects a.*-: 
well as the effects of surface curvature has been neglected 
in the present fom of the turbulence model. Modifications 
to include these additional effects continues to be a subject 
of recent research by different schools. 

Equations (4), (5), (7), (8) and (9) are the equations to 

be solved for w, <l > , H, E and e to study the turbomachine flow 
problem. One can observe that these equations constitute a 
system of coupled elliptic partial differential eauations, 
involving second ordar derivatives. From the nature of the 
problem, none of the terms are negligible in the governing 
equations. In order to solve these elliptic eauations, it is 
necessary to define a selected region in the physical domain 
with boundary conditions specified for all the flow variables 
along the boundaries. The computatio \1 domain used is shown 
in Fig. 2. 

Boundar y Conditions 

Inlet B -jndarv AN : 

Irrotational and prescribed uniform flow conditions are 
imposed along this boundary. Therefore, 


It 

3m 


1 3 ^_ 
r 3® 


tane inlet ■ 


31 . m _1 

34 * 2t/z 


U) 


= - [2£ sina] 


inlet 


H = H . . . 
inlet 


( 10 ) 


The value of the inelt kinatic energy of turbulence, E. , 

inlet 

is specified to be 0.05% of the square of inlet velocity. 

The dissipa k on rate is deduced fr or. the turbulence length 
scale [9] : 


e 


E 


3/2 

inlet 



(10a) 


where 1 is taken to be 0.006 m for the case studied in the present 
work. 
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The Periodic Boundaries AB, NM and FG , IH : 

The periodicity condition is enforced along these boundaries. 
This implies that the value of all fluid properties as well as 
their <j> derivatives, except for the stream function p , have 
equal values at every two corresponding points along AB, NM. 
Regarding the stream function, the value of p's are made to 
differ on the two boundaries AB, NM by a constant value that is 
proportional to the mass flow rate, M. Similarly, the same 
conditions are imposed along FG, IH. 

The Blade Surfaces MI and BF : 

A special approach is needed to compute the flow variables 

near the blade surfaces for the following reasons. In the central 

region of the flow, the gradients of the flow properties are 

usually not very steep, and a moderately fine finite difference 

grid yields accurate solutions. However, close to the blade 

surfaces, the variations of flow properties are much steeper, 

thus requiring an extremely fine grid for accurate computation. 

This implies that the total storage requirement in the computer 

s till be rather large. Moreover, the present form of the 

turbulence model as given by equations (8) and (9) is valid only 

for fully turbulent flows. Modifications are required to make 

it applicable near laminar sublayer regions where the turbulence 

2 

Reynolds number (- pE /y^e) is low. T_ overcome these diffi- 
culties, an approach based on the use of wall functions [9] will 
be employed in the present study. In this approach, all the 
finite-difference gird nodes (except for those representing the 
blade surface) are located in the fully turbulent region. It is 
then assumed that a logarithmic velocity profile exists in the 
region between the blade surface and the point (w+1) located in 
the turbulent region. This can be expressed as follows: 


W 


w+1 


M 


d’p 

W 


= c 


w+1 


!/ 4 e 1/2 
D w+1 


K in[e 


AN p C 


1/4 

D 


gl/2/ i 

V 171 z 1 


(ID 


where the subscript (w+1) is used for the values at the internal 
grid node (w+1), C Q is the constant in the turbulence model, 
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K and e are the log-wall constants (K = 0.42 and e = 9.0), 
and IN is the distance of (w+1) from the blade surface. The 
velocity W w+ ^ is assumed to be parallel to the wall-shear 
stress which is given by: 


1 1 2 

x = p C!' ^ E 
w D w+1 


(11a) 


The value of the energy-dissipation rate at the internal 
grid node is determined from the requirement that the 

length scale varies linearly with the distance from the 
blade surfaces 


-3/4 _3/2 
L D w+1 


(11b) 


The value of the energy of turbulence E w+ -i is calculated from 
the governing equation (8), with the following modifications. 
The diffusion of energy to the blade surface is set equal to 
zero, and the generation term D in equation (8) is modified 
to account for the value of the wall shear stress as given by 
equation (11a) . The dissipation term is also modified 
according to equation (lib) , and is assigned an average value 
for the near-wall node as follows: 


N w+1 

oe = o / edN 
o 


P C 


3/4 

D 



w+1 

/ 



(11c) 


For the vorticity boundary condition, equation (11) is 
used in conjunction with the stream function equation (4) to 
determine w at the grid node (w+1) . 


Downstream Boundary GH : 

The conditions of zero gradients in the meridional direction 
are imposed along this boundary for H, E and e, i.e. 


5H 

^m 


0 


and 



( 12 ) 


■m 


0 


11 


(12a) 



For the vorticity u>, the absolute value is taken to be zero, 
thus 


u> 


[23 


sina] 


exit 


(12b) 


For the stream function y , the downstream flow velocities, 
which may be used to determine o derivatives along GH , and that 
guarantee a unique solution to the problem are unknow apricri. 
Therefore, a supplementary condition based on the conservation 
of angular momentum principle is employed in the present work 
to yield the required unique solution. The details of the 
procedure can be found in reference [7] . 


Coordinate Transformation 

The solution of equations (4) -(9) subjected to the appro- 
priate boundary conditions discussed earlier are carried out 
numerically. In order to reduce the complexity of handling the 
numerics near the curved boundaries of the blade surface, a 
coordinate transformation of the (n,b,n) system to a contracted 
boundary fitted coordinate is required. The overall effect of 
this transformation is to produce a square field in which the 
blade surfaces become straight and parallel. As pointed cut 
earlier in reference [7] , two transformations are used to 
implement the generation of the boundary fitted coordinates for 
the blade -to-blade domain in Fig. 3a. The first one is obtained 
by defining a stretched meridional coordinate, x, given by: 



The second transformation [6] generates the boundary fitted 
coordinates ■' and n through the numerical solution cf the 
following equations for x(^,") and i(F,n). 



23 


3-r-n 


+ y 



Q C , r ) 


5 






PC ,"■) 


(13b) 


1 O 



where 
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(= 21 , 2 + , 11 ) 2 

V 3V ' 3r 


:X 3x 3<t 3(j) 

3n 3 Z 3£ 3n 


= (— ) 
l 3£' 


(Ii\ 


(13c) 


The functions Q and P in the above equations are appropriately 
chosen to provide control over the soacing of the coordinate 
lines in the field. 

Mapping the region of interest in Fig. 3a in terms of the 
new boundary-fitted system of coordinates (S,n) yields a fixed 

J 

square field in the final transformed domain as shown in 
Fig. 3b. Thus, it is possible to carry out all flow calculations 
on the fixed square field, using uniform grid with no inter- 
polation required regardless of the blades shape in the physical 
space . 

Before transforming the flow equations (4) -(9) from the 
physical domain to the fixed square field, it is found con- 
venient to express the flow variables in these ecruations in 
the following dimensionless form: 


uir. 


to = 


( W m ) _ 
m o 


H* = - H 


(W 2 ) 
m o 


* r * b 

r = — , b = — 


E = 


(w ;> 

m o 


* 

£ = 


er . 


P r i 




( wjV 

m o 


Re = (W ) 


m o u 


R = 
o 


(14) 


where the subscripts o, t denotes the conditions at the upstream 
boundaries (A-N) and (M-B) respectively. 

By applying the chain rule for partial differentiation 
and using the expressions given by the relation (14) , the 
governing flow equations in the transformed domain can be written 
in terms of the new variables (tj,n) as [7]: 
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Stream Function Ectuation: 
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Vorticity Transport Equal i on : 
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E nergy Equation : 
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Turbulent Kinetic Energy Equation : 
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Dissipation Rate Equation : 
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★ * 

where z is the number of blades, D and G. are the nondimens ional 
equivalent of the source terms D ana G 0 in equations (7b) and (7a) . 
The transformation parameters, J 6 , {?, y in the above equations 
are defined, in the present work, using Thompson's numerical 
package [6] . 

As pointed out in reference [7], the solution of equations 
(14) -(18) in the (5/H) domain provides the required distribution 
of the flow variables t , oj , H, E and e in the physical space. 

In order to solve these equations numerically, a suitable 
finite difference method is required. The computational method 
used is outlined in the following section. 
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SOLUTION PROCEDURES 


The derivation of the general finite difference equations 
corresponding to equations (14) -(18) is described first. This 
is followed by the procedures employed for the simultaneous 
solution of the resulting difference equations. 

Referring to Fig. 3b, the finite difference representation 
of equations (14) -(18) is obtained using standard five point 
difference operator. Thus, 


'p ■ . = A’p . . + A 0 4 1 . . . + A ’i . • , + A.iJj. . , 

1,1 1 i+l, j 2 i-l, j 3 'i,i+l 4 1,3-1 


+ a 5 , 


(19) 
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+ 


'5 ' 


( 21 ) 


* 

E . 
1/1 


D.E. , . 

1 1 + 1,1 


+ D„E. - . 

1 1 - 1,1 


+ D-, E ■ . i +D.E. . . 

3 1,1+1 4 1,1-1 


and 


k k k k k 

£ . • = K . £ ■ , . . + K~£ • . • + K, £ . + K . S . . . 

1,1 1 1+1,1 2 l-l/l 3 1,1+1 4 1,1-1 


( 22 ) 


K r 


(23) 


where the value of the different coefficients in the above 
equations are given in Table 1. A combination of central and 
upwind differences are used for evaluating these coefficients. 
The details of the differencing techniques is reported in 
reference [ 7 ] . 
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Solution of equations (19)- (23) will be obtained using a 
successive relaxation iterative method. In using this method 
there are alternate methods for handling the coupling between 
the governing equations. For example, if these equations are 
solved in sequence, the coupling has to be restored through 
repeated iterations. Not all iterative schemes, however 
converge to a solution and experience proves that the solution 
of a system of coupled nonlinear equations converges for 
laminar flows, but diverges for turbulent flows. The problem 
of instability usually arises from the coupling between the 
mean velocity profile and effective viscosity profile. 

Reference [12] suggests the use of a successive convergence 
technique to retard and dampen this coupling to ensure con- 
vergence under any set of boundary conditions. 

The successive convergence technique is implemented in 
the present work by decoupling the turbulent kinetic energy 
and dissipation rate equations (22) and (23) from the remaining 
equations. The former are solved sequentially rather than 
simultaneously with the vorticity and the stream function 
equations (19) and (20). The details of the procedure is as 
follows. The stream function and vorticity equations are 
first solved with a frozen viscosity field using the same 
numerical scheme outlined previously in reference [7] . Once 
a convergence solution has been obtained, the kinetic energy of 
turbulence and the dissipation rate equations are then solved 
using a Gauss-Siedel method to obtain a new round of frozen 
viscosity field which will be used in the next cycie of 
convergent solution for p, o . The step described immediately 
above is executed several times until a convergence criteria 
is satisfied. When the calculations have reached this stage, 
the energy equation (21) is then solved using again the Gauss- 
Siedel method. 

The computer program used to carry out the actual flow 
calculations is an expanded version of the one used in Part I 
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of this study. In order to keep the computer time within 
reasonable limits (usually less than nine minutes for turbulent 
cases on an AMDAHL 470) , the unit square region has been 
divided into 30 step sizes in n direction and 40 in the t, 
direction. These grids have not been optimized, however 
they will suffice for the present purpose, which is to 
demonstrate the capability of carrying out calculations of this 
type of flow. 


3. COMMENTS ON THE COMPUTATIONAL METHOD 

Before proceeding to a detailed discussion of the numerical 
solution of the flow case studied, a few comments are given 
first. 

Special consideration had to be taken into account during 
the numerical generation of the boundary fitted coordinate 
system for the blade- to-biaae domain of Fig. 3. This pertained 
to the enforcement of the periodic boundary conditions in an 
effective manner. In order not to resort to additional 
iterations solely for the purpose of enforcing the periodicity 
conditions, it was required that the t, coordinate line emanating 
from one point on the lower boundary AB of Fig. 3 must meet 
the upper boundary NM at the point having the same meridional 
location. Generally, this constraint was easily satisfied, 
since as explained in reference [7], the assignment of the 
^-values to the boundary points via the functions £. , £ g^ 
and g 2 is arbitrary. Although the distribution in the Z lines, 
in these cases, did not necessarily produce orthogonal coor- 
dinates, as exemplified in Ficj. 3, the lack of orthogonality 
created no problem at all. 

Some problems arose with the convergence characteristics 
of the kinetic energy of turbulence equation during the numerical 
solutions. Guided by an analysis, for a simple test case with 
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few grid points, and by subsequent experience on the computer, 
the cause of the problem was traced and defined to be the source 
term (D) in equation (8). Since this term is a function of 
the square of the kinetic energy itself [see equations (7b) 
and (7c)], it might therefore be anticipated that this term 
vary so widely, during the iteration procedure, so as to 
provoke divergence. To reduce these variations, a rearrangement 
of the finite difference equation (22) is sought. For this 
the source term, D, was written as 


D = ^e F l = ^ F 1 

+ E (C D 

P x- 

F l } 


(24a) 

where 






3W 2 

P i- J ' ( 5r> + 

. 3W . 
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r 

sinct) 
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3W, 

+ ( ' r 
am 

, aw 

1 m 

r 

W A 

r 

sina) 

2 

(24b) 


Thus, by substituting equation (24a) in equation (22) , and 
rearranging, one obtains: 


* D. E . . . +D_E. , . + D,E . ... + D ,E . . . + D c 

__ 1 i + l,j 2 i l_,j 3 1,^+1 4 _i,jgl 5 

lb * 




(1 . ac, an C F* E*/c*/D). ■ 

u i 01,3 


(24c) 


where 


D_ = Li an [ -y ,F. 
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2’T 


U C 


*2 
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0 t 
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he s b n 

Llj 


* 

-r~) 


|n ^Re S H ^ 

lb 


( 24d) 


With this new substitution formula, the variations in the 
modified source term was found to be much loss than those t 
the original one and convergence was restored for all case; 
studied in the present work. 



A problem similar to that immediately discussed above was 
developed with the source term in the dissipation rate equation 
(18). Again, the finite difference equation (23) was rearranged 
in a similar manner so as to dampen the variation of the source 
during the numerical solution. 

It is appropriate to mention at this stage a special 
practice was employed near the blade exit region along the 
boundary FI (see Fig. 2). Under certain circumstances , a 
smooth variation for the kinetic energy of turbulence and its 
dissipation rate, was not achievable near the blade exit. 

This in turn caused some difficulties in the calculations, 
ana eventually induced divergence especially in cases involving 
small convergence criteria. In order to alleviate this kind of 
restriction, the condition given by equation (12a) was applied 
along the boundary FI in Fig. 2 instead of the boundary HG. 
Although this modification ensured convergence, under any set 
of conditions, it is felt that further study is recommended 
to remove the inadequacy resulting from the neglect of the decay 
of E and e downstream the blade row. 

The calculation procedures used to establish convergence 
for each cycle of solution for the stream function and vorticity 
were covered in some detail in Volume I of this study [7]. 
Typically, intermediate convergent solutions at given effective 
viscosity distributions were obtained after 50-30 iterations on 
the average. A succession of 15-25 rounds of intermediate con- 


vergent solutions with frozen effective viscosity fields were 
required for a complete numerical solution of a typical problem. 
In the case studied hero, the determining convergence criteria, 
in each cycle of w and i calculations, turned out to be the 
vorticity tolerance. That is, when the vorticity tolerance was 
finally satisfied, the stream function error was considerably 
below the required values. Likewise, the determining convergence 
criteria during the subcycle or. E and •; calculations turned 
out to be the kinetic enure v of turbulence tolerance . 


2 U 



4. RESULTS AND DISCUSSION 


The procedure described in the previous sections is employed 
to obtain the flow details in a cooled radial inflow turbine 
rotor. The rotor geometry is shown in Fig. 4 together with the 
shape and thickness distribution of the computational surface S^. 
The blade-to-blade shape in the physical domain as well as the 
boundary fitted coordinates employed in the solution are those 
shown in Fig. 3. Additional summary data related to the blades 
geometry and the configuration of the computational surface S^ 
can be found in Table 2. The following operating conditions for 
the turbine are used in the analysis: 

Turbine inlet total temperature, T*. , 1083 °K 

z 3 

Turbine inlet density, p, 1.0060 Kg/m 

Rotational speed, Q, 38,500 r.p.m. 

Rotor inlet flow angle, 8^ n , 62.5° 

Rotor exit flow angle, 8 ex ^ t » 63° 

Meridional component of the relative velocity at rotor 

inlet, W m , 66.2 m/s pW r 

Flow Revnolds number at rotor inlet, Re (= — • . ., 

3.34 x 105 “e lnlet 

Prandtl number, 0.8 

The flow through the turbine is considered to be incompressible 
and the blade surfaces are cooled and kept at a uniform tem- 
perature of 550 e K. 

Results of the analysis are presented in Figs. 5 -through 9 
as contour plots for the distribution of the stream function, 
the velocity, the kinetic energy of turbulence as well as the 
temperature within the blade passages. 

The stream function contour plots are givr*n in Fig. 5. 

The flow contours depicted reflect the appearance of a recir- 
culating eddy near the pressure surface of the blades. The flow 
rate within the eddy amounts to about 2.5t of the flow rate 
through the turbine (as indicated by the designation of the 
streamlines). This low flow rate implies low relative velocities 
within the eddy as will be discussed shortly. The manner in 
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which the streamline (v ■ 1) leaves the blade surface indicates 
that a separation bubble is formed over the aft portion of the 
blade suction surface. This behavior seemed to be caused by 
the inability of the flow to accommodate the large blade deflection 
angle in this region. 

The distribution of relative velocity, W, across the 
blade-to-blade passage, nondimens ionalized with respect to the 
blade tip speed, is shown in Pig. 6, The contours show 

a core of high relative velocity centered at a point located 
at approximately 15 percent of the passage width from the blade 
suction surface. A region of distinct velocity deficiency is 
also observed in the upper half of the passage extending to 
about 70 percent the blade chord and is most pronounced near 
the blade pressure surface. This region coincides with the 
location of the recirculatory eddy shown in Fig. 5. The data 
in Fig. 6 also indicates that the flow experiences a high 
acceleration rate near the blade pressure surface as it approaches 
the passage exit. This behavior is similar to that of flow 
approaching a turning duct and is expected. 

In order to gain some insight into the development of 
the boundary layer over the blade surfaces, the velocity profiles, 
across the rotor passage, at different meridional locations 
are introduced in Fig. 7. The locations are selected to 
correspond to the meridional distances 1.0, 4.0 and 6.0 cms 
respectively. The profiles of Fig. 7a and 7b show reversed flow 
regions existing near the blade pressure surface and steep 
gradients in the velocity near the suction surface. There it a 
noticeable change in the velocity profiles at subsequent do. '-.stream 
stations. A jet type of profile exists near the turbine exit as 
illustrated in Fig. 7c. The complex flow pattern observed 
generally in Figs. 6 and 7 indicates the significant extent 
of the complex interaction mechanism that occurs between the 
different flow regimes (i.e. boundary layers, recirculating 
eddies and separation regions, etc.) within the rotor channel. 

As a consequency of such interaction, ru clear distinction can be 
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made between the boundary layer regions and the inviscid core. 

This result emphasize!* the inaccuracies involved in using a 
nonviscous flow solution, such as the one shown in Fig. 7 (with 
dotted lines) , to generate boundary layer characteristic 
parameters for flows in turbomachine rotors. 

Figure 8 shews the distribution oi the kinetic energy, 

of turbulence between the blades. The turbulence intensities 

« 2 2 

are piottad using the nondimensional quantity E » E/f. r fc ^ . 

Near the passage inlat and up to 60 percent of the blade chord 
peak turbulence intensities are observed to occur in the mixing 
region of the recirculating eddy noted in Fig. 5. The flow 
energy dissipated in this process of recirculation is responsible 
for the velocity defficiency encountered near the blade pressure 
surface as shown in Fig. 5. On the whole, it could be recognized 
that the turbulence intensities are generally higher than those 
encountered in a stationary passage. The maximum intensities 
are seen to occur near the blade surfaces and especially at the 
passage exit. The relatively high values of kinetic energy 
of turbulence shown in the figure away from the blade surfaces 
reflect the extent of turbulent mixing even near the mid passage. 
The classical assumption that the viscous and turbulence effects 
are confined to very thin regions near the blade surfaces is 
evidently inaccurate when dealing with flow in turbomachinery. 

It is with the present type of analysis that involves the solution 
of the full field equations in addition to a suitable closure 
model for turbulence that a realistic flew prediction can be made. 
It should be emphasized, however, that the kinetic energy of 
turbulence results presented here must be viewed as qualitative. 
This is attributable to the fact that the additional production 
of turbulence energy by the rotational effects as well as the 
effects of blades curvature has been neglected in the present 
form of turbulence model. 

Finally, the static temperature distribution within the 
blade-to-blade passage, nondimensional ized with respect to 
the inlet total temperature is given in Fig. 9. Inspection of 
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the contours plot indicates that the steepest gradients in 
temperature profiles occur near the blade suction surface at 
the channel inlet and near the pressure surface at the channel 
exit. This behavior is plausible, since the velocities in these 
regions are high (see Fig. 6) . Conversely, it is observed that 
the thermal layers tend to thicken on the suction surface at 
trailing edge and also near the pressure surface where the 
recirculating zone exists. Since the rate of heat transfer 
to the blades is proportional to the temperature gradients near 
the surface, therefore, one can conclude that maximum rate of 
heat transfer occurs at the blade leading edge and also near 
the blade suction surface. This result has important impact 
on the process of designing the cooling passages for new 
impellers, particularly when the design procedure is based on 
specifying the blade temperature distribution. The present 
results suggest that in order to keep the blades cooled at 
constant temperature, it is necessary to place large number of 
cooling passages near the blade leading edge and also along 
the suction surface in the neighborhood of the impeller inlet. 

5. CONCLUSIONS 

A method has been developed which is suitable for analyzing 
the viscous flow within turbomachine passages. The field 
analysis was based upon the numerical integration of the full 
incompressible Navier-Stokes equations over a predetermined 
set of blade-to-blade stream surfaces. The method employs a 
nonorthogonal body-fitted coordinate system to account for the 
most complicated blade geometries. Effects of turbulence were 
modelled using a two-equation turbulence model. 

Significant viscous flow results were obtained by applying 
the method of analysis to a radial inflow turbine. The 
predicted flow patterns indicated the severity of the complex 
interaction mechanism that occurs between the different flow 
regimes (i.e. boundary layers, recirculating eddies, separation 
zones, etc.). This result emphasized clearly the inaccuracies 
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involved in using a nonviscous flow solution with classical 
boundary layer techniques to obtain viscous flow details 
within turbomachine rotors. Heat transfer results pointed 
out that in order to keep the rotor blades cooled at uniform 
temperatures, it was necessary to place a large number of 
cooling passages near the blades leading edge and also along 
the suction surface in the vicinity of rotor inlet. 

On the whole, it is recognized from the results presented 
here and in Volume I, that the present flow analysis provides 
a good prediction of the actual flow behavior within the passage 
of turbomachine rotors. This accomplishment cannot be 
attributed to any single element of the analysis, but results 
from the combinations of many factors. The major one is believed 
to be the preservation of the ellipticity of the problem by 
working with the full Navier -Stokes equations. Other important 
factors include the transformation, discretization and the 
numerical procedures used that provided stable calculations 
of good accuracy. 

As with all research efforts of this magnitude, several 
areas for future work are recognized during the course of the 
present study. These include the use of variable grid in the 
numerical scheme, the inclusion of the effects of rotation and 
surface curvature in the turbulence model and the inclusion of 
the flow compressibility in the solution procedure. Similarly, 
the inaccuracies involved in using the present method to 
analyze flows including three dimensional separated zones are 
identified and documented (see Vol. I). 
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NOMENCLATURE 


A 

B 


1 

1 


b 

C 

C 


1 

1 





D 

D 


1 


E 

G 

h 




» 


H 

J 

K r K 2 ,K 5 

l 

k 


n 

AN 

P 

Pr 

r 

R 

Re 

S CE 

S C£ 


T 

W 

W 

W 


m 


W, 


Coefficients of the finite difference equations. 
Coefficients of the finite difference equations. 

Normal stream annulus thickness, m. 

Coefficients of the finite difference equations. 
Constants in the turbulence model. 

Specific heat at constant pressure, J/(Kg) (°K) . 
Dissipation function. 

Coefficients of the finite difference equations. 
Kinetic energy of turbulence, J/Kg. 

Denoting source terms in the flow governing equations. 
Static enthalpy, J/Kg. 

Total enthalpy, J/Kg. 

Jacobian matrix. 

Coefficients of the finite difference equations. 
Mixing length, ra. 

Mass flow per blade flowing through the stream 
annulus, Kc-/sec. 

Outward unit normal to the stream surface, S^. 

Distance of the near-wall grid point (w+1) from the 
blade surface. 

2 

Static pressure, N/m or blade pitch. 

Effective turbulent Prandtl number. 

Radius from axis of rotation, m. 

Universal gas constant, U/(Kg) (°K) . 

Reynolds number. 

Schmidt number for kinetic energy of turbulence. 
Schmidt number for dissipation of kinetic energy of 
turbulence . 

Temperature, °K. 

Relative velocity vector, m/sec. 

Magnitude of W, m/sec. 

Meridional component of the relative velocity vector, 
m/sec. 

Tangential component of the relative velocity vector, 
m/sec. 
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X 


Stretched meridional coordinate, 
z Axial coordinate/ m. 

Z Number of blades. 

a Angle between m and z, rad., see Fig. 1. 

6 Coordinate transformation parameter, or angle between 

relative velocity vector and meridional plane, rad., 
see Fig. 2. 

Y Coordinate transformation parameter. 

6 Coordinate transformation parameter, 

p Fluid density, Kg/m^ . 

2 

y Effective viscosity, m /sec. 

e 2 
y^ Laminar viscosity, m /sec. 

ft Rotational speed, rad/sec. 

4) Relative angular coordinate, rad. 

\p Stream function. 

oj Vorticity, 1/sec. 

e Dissipation of kinetic energy of turbulence, J/Kg. 

C Boundary fitted coordinate, 

n Boundary fitted coordinate. 

Superscripts 

Mean value. 

* Denotes nondimens ional quantity. 


Subscripts 

e 

i 

m 

w 

w+1 

if j 
inlet 
exit 
total 
t, tip 


Effective 

Laminar. 

Meridional component 
Wall value. 

Pertaining to points in the flow field at a distance 
An away from w. 

Tangential component. 

Denotes field position in (£,n) domain, see Fig. 3b. 
Inlet or upstream. 

Exit or downstream. 

Total conditions. 

Rotor tip. 
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TABLE 1 - COEFFICIENTS OF 


THE FINITE DIFFERENCE EQUATIONS 


t 


A 1 = 26 An/AC/A 0 /(b* +1 + b*) 
A 2 = 26An/A£/A o /(b*_ 1 + b*) 


A 3 = y A£/An/A Q /b* 

A 4 = Y A£/An/A Q /b* 

(K |f) i/a 0 

A = A. + A- + A, + A. 
o 1 2 3 4 


r t * 2 2 

A 5 = A£ An [~ U) R J z/2tt 
o o 


e i-r— 
- H ( b * 


ii* <L_ 
3n' - 3n 


where 


t All unsubscript quantities in this table are evaluated at the grid 
node (i,j). Note that the stream channel thickness b* is constant 
for all j. 1 



TABLE 1 - CONTINUED 


B 1 “ l B ll + B 12 b l«'< 5 e >i+l,j 1 b X + l ' 


B 2 * ' B 21 + B 22 V b l-1 


B 3 tB 31 + B 32 b i /<Be) i, j+1 1 b i/ b i+l 


B. = B., + B., b./(Re. . ,) 

4 41 42 x 1 , 3-1 


where 


B 11 

- J< *n 

+ 1*11 

B 12 

= 6 An/A?/B Q 

B 21 

* J <*22 

+ 1*22 

B 22 

= 6 An/A£/B Q 

B 31 

‘ J( *33 

+ 1*33 


)/8B 0 

)/8B 0 


)/8B o 


B32 = Y * 5 /An/B o 
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®41 = J( *44 + I <*44 I 


B 42 = Y AC/A n/B 
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2 * R 

'j * r a*, m ^ w \ j t o 3 /Sina 3ty % , 3 / s i not d ^ \ i /n 

2 AC An b [0 (^-) + J — — (— *- ^ (— *- 


Re 


Re 3C ' 3n 3n ' b * 3C "o 


* ~ 


B o = B 11 + B 21 + B 31 + B 41 + (B 12 + B 22 + B 32 + B 42> b / Re 


*11 " *i+l,j-l + *i, j-1 " *i+l , j+1 " *i,j+l 


*22 “ *i 


-l,j+l + *i , j+1 “ *i-l,j-l ” *i , j-1 


*33 *i+l 


,j+l + *i+l , j “ *i-l , j+1 “ *i-l , j 


*44 *i-l,j-l + *i-l,j ~ *i+l,j-l - *i+l,j+l 
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(C 


21 


(C 


31 


(C 


41 


+ C 


12 


+ C 


22 


+ C 


32 


+ C 


42 


)/C o 

)/C o 

)/C o 

>/ C o 


where 

C 11 = J b o ( ^11 + Un l)/ 0b * 

C 12 = 6 An/A£/Pr/[ (Re) i+1> j + <Re) ± ^ ) /2 . 0 


C 21 = J b o <*22 + l^2 2 l )/8b * 

C o0 = 6 An/AC/Pr/[ (Re) + (Re). .]/2.0 

Z> £ i-L f j X / J 


* The Prandtl Number is assumed to be constant in this table 
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TABLE 2: BLADE GEOMETRY DATA RADIUS 


Blade Lower Surface Data 


Meridional 
Distance m(m) 

Angular 

Coordinate 

0 ! 62539D-03 

0.56273 

0 . 86000D-G2 

0.55900 

0.16000D-01 

0.55560 

0 . 23500D-01 

0.55220 

0.29030D-01 

0.55050 

0 . 34280D-01 

0.54900 

0 . 39540D-01 

0.54640 

0 . 46230D-01 

0.51540 

0 . 5394CD-01 

0.42600 

0 .6197 0D-01 

0.23650 

0 . 67011D-01 

0.23898D-01 


Blade Upper Surface Data 


Meridional 
Distance m(m) 

Angular 

Coordinate 

4> 

0 . 62539D-03 

0 . 84710D-02 

0 . 86000D-02 

0.12200D-01 

0 . 16000C-01 

0 . 15600D-01 

0 . 23500D-01 

0 . 18900D-01 

0.29030D-01 

0 . 20900D-01 

0 . 34280D-01 

0.21500D-01 

0 . 39540D-01 

0.16000D-01 

0 . 46230D-01 

-0 . 11600D-01 

0 . 53940D-01 

-0 . 96200D-01 

0.61970D-01 

-0.27910 

0 . 68363D-01 

-0.53070 


Derivative 

2nd 

Derivative 

-0.45678 

-7.2554 

-0.46035 

6.3594 

-0.43120 

-11.935 

-0.35261 

45.238 

-0.34501 

•4 .542 

-0.52816D-01 

154.86 

-1.7845 

-813.31 

-7.6118 

-928.77 

-15.982 

-1242.5 

-33.844 

-3206.4 

-50.743 

-3497.8 


Derivative 

2nd 

Derivative 

0.45678 

5.7000 

0.46655 

-3.2487 

0.45729 

0.74568 

0.40262 

-15.325 

0.32213 

-13.787 

-0.26520 

-209.96 

-2.0543 

-470.30 

-8.6949 

-917.03 

-15.993 

-1495.0 

-30.342 

-0278.9 

-50.743 

-4303.7 
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TABLE 2: BLADE GEOMETRY DATA 


Meridional 
listance m(m) 

Radius from 
Center of 
Rotation 

-0.007079 

0.083527 

-0.003684 

0.080123 

0.0 

0.076450 

0.003498 

0.072986 

0.006853 

0.069697 

0.010048 

0.066603 

0.013104 

0.063701 

0.016034 

0.060999 

0.018845 

0.058511 

0.021535 

0.056258 

0.024127 

0.054236 

0.026630 

0.052443 

0.029055 

0.050874 

0.031413 

0.049519 

0.033711 

0.048348 

0.035959 

0.047338 

0.038162 

0.046477 

0.040329 

0.045762 

0.042466 

0.045177 

0.044574 

0.044690 

0.046662 

0.044266 

0.048730 

0.043881 

0.050782 

0.043516 

0.052817 

0.043155 

0.054835 

0.042782 

0.056835 

0.042392 

0.058816 

0.041987 

0.060776 

0.041570 

0.062717 

0.041145 

0.064639 

0.040729 

0.066544 

0.040353 

0.068434 

0.040051 

0.070308 

0.039843 

0.072176 

0.039709 

0.074039 

0.039618 

0.075398 

0.039550 

0.077754 

0.039499 

0.079607 

0.039462 


sinct 

Stream Annulus 
Thickness, b(m) 

-1.004280 

0.000967 

-1.000120 

0.000956 

-0.993984 

0.000960 

-0.985732 

0.000982 

-0.974760 

0.001013 

-0.960602 

0.001043 

-0.937357 

0.001068 

-0.904973 

0.001090 

-0.862745 

0.001110 

-0.810406 

0.001128 

-0.749061 

0.001144 

-0.681938 

0.001156 

-0.610477 

0.001160 

-0.540900 

0.001155 

-0.478504 

0.001157 

-0.420336 

0.001179 

-0.360370 

0.001216 

-0.299578 

0.001253 

-0.250202 

0.001234 

-0.214878 

0.001310 

-0.192922 

0.001335 

-0.180691 

0.001362 

-0.176406 

0.001389 

-0.179872 

0.001416 

-0.190020 

0.001441 

-0.200095 

0.001465 

-0.208736 

0.001490 

-0.215980 

0.001514 

-0.220585 

0.001540 

-0.209999 

0.001568 

-0.131434 

0.001595 

-0.135325 

0.001620 

-0.088592 

0.001642 

-0.057745 

0.001660 

-0.042105 

0.001673 

-0.331892 

0.001681 

-0.023363 

0.001686 

-0.016506 

0.001688 
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FIG. 5b. GRID ORDERING SYSTEM IN THE TRANSFORMED DOMAIN. 
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. f|, HUB-SHROUD PROFILE WITH STREAM SURFACE Sp USED FOR 
BLADE-TO-BLADE ANALYSIS. 
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flOMDIMEMS TONAL VELOCITY DISTRIBUTION AT DIFFERENT 
MERIDIONAL LOCATIONS. 
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